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Abstract A remarkable number of different numerical algorithms can be understood and 
analyzed using the concepts of symmetric spaces and Lie triple systems, which are well 
known in differential geometry from the study of spaces of constant curvature and their 
tangents. This theory can be used to unify a range of different topics, such as polar-type 
matrix decompositions, splitting methods for computation of the matrix exponential, com- 
position of selfadjoint numerical integrators and dynamical systems with symmetries and 
reversing symmetries. The thread of this paper is the following: involutive automorphisms 
on groups induce a factorization at a group level, and a splitting at the algebra level. In this 
paper we will give an introduction to the mathematical theory behind these constructions, 
and review recent results. Furthermore, we present a new Yoshida-like technique, for self- 
adjoint numerical schemes, that allows to increase the order of preservation of symmetries 
by two units. Since all the time-steps are positive, the technique is particularly suited to stiff 
problems, where a negative time-step can cause instabilities. 
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1 Introduction 

In numerical analysis there exist numerous examples of objects forming a group, i.e. ob- 
jects that compose in an associative manner, have an inverse and identity element. Typical 
examples are the group of orthogonal matrices or the group of Runge-Kutta methods. Semi- 
groups, sets of objects close under composition but not inversion, like for instance the set of 
all matrices and explicit Runge-Kutta methods, are also well studied in literature. 

However, there are important examples of objects that are neither a group nor a semi- 
group. One important case is the class of objects closed under a 'sandwich-type' product, 
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(a, b) i — y aba. For example, the collection of all symmetric positive definite matrices and all 
selfadjoint Runge-Kutta methods. The sandwich-type composition aba for numerical inte- 
grators was studied at length in 117] and references therein. However, if inverses are well 
defined, we may wish to replace the sandwich product with the algebraically nicer symmet- 
ric product (a.b) H- ab~ l a. Spaces closed under such products are called symmetric spaces 
and are the objects of study in this paper. There is a parallel between the theory of Lie groups 
and that of symmetric spaces. For Lie groups, fundamental tools are the Lie algebra (tangent 
space at the identity, closed under commutation) and the exponential map from the Lie al- 
gebra to the Lie group. In the theory of symmetric spaces there is a similar notion of tangent 
space. The resulting object is called a Lie triple system (LTS), and is closed under double 
commutators, [X, [Y,Z]]. Also in this case, the exponential maps the LTS into the symmetric 
space. 

An important decomposition theorem is associated with symmetric spaces and Lie triple 
systems: Lie algebras can be decomposed into a direct sum of a LTS and a subalgebra. The 
well known splitting of a matrix as a sum of a symmetric and a skew symmetric matrix is an 
example of such a decomposition, the skew symmetric matrices are closed under commuta- 
tors, while the symmetric matrices are closed under double commutators. Similarly, at the 
group level, there are decompositions of Lie groups into a product of a symmetric space and 
a Lie subgroup. The matrix polar decomposition, where a matrix is written as the product of 
a symmetric positive definite matrix and an orthogonal matrix is one example. 

In this paper, we are concerned with the application of such structures to the numeri- 
cal analysis of differential equations of evolution. The paper is organised as follows: in §2 
we review some general theory of symmetric spaces and Lie triple systems. Applications 
of this theory in numerical analysis of differential equations are discussed in §3, which, 
in turn, can be divided into two parts. In the first (§3. 1 — §3.3), we review and discuss the 
case of differential equations on matrix groups. The properties of these decompositions and 
numerical algorithms based on them are studied in a series of papers. Polar-type decompo- 
sitions are studied in detail in 1211 . with special emphasis on optimal approximation results. 
The paper [32] is concerned with important recurrence relations for polar-type decompo- 
sitions, similar to the Baker-Campbell-Hausdorff formula for Lie groups, while |34, 10, 12] 
employ this theory to reduce the implementation costs of numerical methods on Lie groups 
[9 20]. We mention that polar- type decompositions are also closely related to the more spe- 
cial root-space decomposition employed in numerical integrators for differential equations 
on Lie groups in [23 1 . In 1 13] it is shown that the generalized polar decompositions can be 
employed in cases where the theory of |23 | cannot be used. 

In the second part of §3 ( p.4| and beyond), we will consider the application of this theory 
to numerical methods for the solution of differential equations with symmetries and revers- 
ing symmetries. By backward error analysis, numerical methods can be thought of as exact 
flows of nearby vector fields. The main goal is then to remove from the numerical method the 
undesired part of the error (either the one destroying the symmetry or the reversing symme- 
try). These error terms in the modified vector field generally consist of complicated deriva- 
tives of the vector field itself and are not explicitly calculated, they are just used formally for 
analysis of the method. In this context, the main tools are compositions at the group level, 
using the flow of the numerical method and the symmetries/reversing symmetries, together 
with their inverses. There is a substantial difference between preserving reversing symme- 
tries and symmetries for a numerical method: the first can be always be attained by a finite 
(2 steps) composition (Scovel projection [27 1), the second requires in general an infinite 
composition. Thus symmetries are generally more difficult to retain than reversing symme- 
tries. For the retention of symmetry, we review the Thue-Morse symmetrization technique 
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for arbitrary methods and present a new Yoshida-like symmetry retention technique for self- 
adjoint methods. The latter has always positive intermediate step sizes and can be of interest 
in the context of stiff problems, which typically require step size restrictions. We illustrate 
the use of these symmetrisation methods by some numerical experiments. 
Finally, Section 4 is devoted to some concluding remarks. 



2 General theory of symmetric spaces and Lie triple systems 

In this section we present some background theory for symmetric spaces and Lie triple sys- 
tems. We expect the reader to be familiar with some basic concepts of differential geometry, 
like manifolds, vector fields, etc. For a more detailed treatment of symmetric spaces we refer 
the reader to |4| and [ 14 1 which also constitute the main reference of the material presented 
in this section. 

We shall also follow (unless otherwise mentioned) the notational convention of |4): in 
particular, M is a set (manifold), the letter G is reserved for groups and Lie groups, gothic 
letters denote Lie algebras and Lie triple systems, latin lowercase letters denote Lie-group 
elements and latin uppercase letters denote Lie-algebra elements. The identity element of a 
group will usually be denoted by e and the identity mapping by id. 



2.1 Symmetric spaces 

Definition 1 (See 1 14]) A symmetric space is a manifold M with a differentiable symmetric 
product ■ obeying the following conditions: 

(i) x ■ x = x, 

(ii) x- (x-y) =y, 

(iii) x- (yz) = (x-y) ■ (x-z), 

and moreover 

(iv) every x has a neighbourhood U such that for all y in U x ■ y = y implies y = x. 

The latter condition is relevant in the case of manifolds M with open set topology (as in the 
case of Lie groups) and can be disregarded for sets M with discrete topology: a discrete set 
M endowed with a multiplication obeying (i)-(iii) will be also called a symmetric space. 

A pointed symmetric space is a pair (M, o) consisting of a symmetric space M and a 
point o called base point. Note that when M is a Lie group, it is usual to set o = e. Moreover, 
if the group is a matrix group with the usual matrix multiplication, it is usual to set e = I 
(identity matrix). 

The left multiplication with an element x 6 M is denoted by S x , 

S x y = x-y, VyeM, (1) 

and is called symmetry around x. Note that S x x = x because of (i), hence x is fixed point 
of S x and it is isolated because of (iv). Furthermore, (ii) and (iii) imply S x is an involutive 
automorphism of M, i.e. = id. 

Symmetric spaces can be constructed in several different ways, the following are impor- 
tant examples: 
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1 . Manifolds with an intrinsically defined symmetric product. As an example, consider the 
n-sphere as the set of unit vectors in The product 

x ■ y = S x y = (2xx T — I)y 

turns this into a symmetric space. The above operation is the reflection of points on a 
sphere. This can be generalized to m-dimensional subspaces in W (m < n) in the fol- 
lowing manner: Assume that x = [xi, . . . ,x m ] is a full rank matrix. Define P x the orthog- 
onal projection operator onto the range of x, P x = x(x*x) x*. Consider the reflection 
R x = 2P X — I. Define x ■ y = R x y. This operation obeys the conditions (i)-(vi) whenever 
x,y,z are m x n full rank matrices. In particular, note that (i) is equivalent to R 2 = I, i.e. 
the reflection is an involutive matrix. 

2. Subsets of a continuous (or discrete) group G that are closed under the composition x ■ 
y = xy~ l x, where xy is the usual multiplication in G. Groups themselves, continuous, as 
in the case of Lie groups, or discrete, are thus particular instances of symmetric spaces. 
As another example, consider the set of all symmetric positive definite matrices as a 
subset of all nonsingular matrices, which forms a symmetric space with the product 

a- b = ab~ a. 

3. Symmetric elements of automorphisms on a group. An automorphism on a group G is a 
map a : G — > G satisfying a(ab) = o(a)o(b). The symmetric elements are defined as 

^ = {geG : o(g)=g- 1 }. 

It is easily verified that s& obeys (i)-(iv) when endowed with the multiplication x-y = 
xy~ l x, hence it is a symmetric space. As an example, symmetric matrices are symmetric 
elements under the matrix automorphism a(a) = cT T . 

4. Homogeneous manifolds. Given a Lie group G and a subgroup H, a homogeneous man- 
ifold M = G/H is the set of all left cosets of H in G. Not every homogeneous manifold 
possesses a product turning it into a symmetric space, however, we will see in Theo- 
rem[T]that any connected symmetric space arises in a natural manner as a homogeneous 
manifold. 

5. Jordan algebras. Let a be a finite-dimensional vector space with a bilinear multiplica- 
tiorQx * y such that 

x *y = y * x, x * (x 2 * y ) = x 2 * (x * y) 

(powers defined in the usual way, x m = x * x m ~ 1 ), with unit element e. Define L x (y ) = x * y 
and set P x = 1l} x — L x 2 . Then, the set of invertible elements of is a symmetric space 
with the product 

x-y = P x (y- 1 ). 

In the context of symmetric matrices, take x*y = ^(xy + yx), where xy denotes the 
usual matrix multiplication. After some algebraic manipulations, one can verify that the 
product x-y = P x {y~ l ) = 2x* (jc*y~') — (x*x) *y~' = xy~ l x as in example^! 

Let G be a connected Lie group and let a be an analytic involutive automorphism, i.e. 
CT ^ id and a 2 = id. Let G a denote fixe = {g E G : o{g) — g}, G° its connected component 
including the base point, in this case the identity element e and finally let K be a closed 
subgroup such that Gf C K C G a . Set G a = {x € G : o(x) = jiT 1 }. 



Typically non-associative. 
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Theorem 1 (| 14 1) The homogeneous space M = G/K is a symmetric space with the product 
xKyK = xa(x)~ l o(y)K and G a is a symmetric space with the product x-y — xy~ l x. The 
space of symmetric elements G a is isomorphic to the homogeneous space G/G a . Moreover, 
every connected symmetric space is of the type G/K and also of the type G a - 

The interesting consequence of the above theorem is that every connected symmetric 
space is also a homogeneous space, which implies a factorization: as coset representatives 
for G/G a one may choose elements of G a , thus any x 6 G can be decomposed in a product 
x = pk, where p 6 G a and k 6 G° . In other words, 

x = pk, o(k)=k, o(p)=p~ l , (group factorization). (2) 

The matrix polar decomposition is a particular example, discussed in j |3.1| 

The automorphism a on G induces an automorphism on the Lie algebra Q and also a 
canonical decomposition of Q. Let Q and t denote the Lie algebras of G and K respectively 
and denote by do~ the differential of o~ at e, 

do(X) = -^\ o-(exp(rX)), VX e 0. (3) 
dt \t=o 

Note that da is an involutive automorphism of Q and has eigenvalues ±1. Moreover, Xgt 
implies da(X) = X. Set p = {X e Q : da{X) = -X}. Then, 

Q = p(Bt, (4) 

(4). It is easily verified that 

[t,p]Cp, [p,p]C*, (5) 

that is, t is a subalgebra of Q while p is an ideal in {?. Given X € Q, its canonical decompo- 
sition petisX = P + K, with P € p and K € t, 

X = P + K, de{K)=K, da(P) = -P, (algebra splitting). (6) 

We have already observed that there is close connection between projection matrices, 
reflections (involutive matrices) and hence symmetric spaces. In a linear algebra context, this 
statement can be formalized as follows. Recall that a matrix TI is a projection if TI 2 = IT. 

Lemma 1 To any projection matrix TI there corresponds an involutive matrix S = I — 2TI. 
Conversely, to any involutive matrix S there correspond two projection matrices TI^ = \(J~ 
S) and TI^ = \ (/ + S). These projections satisfy TI^ + TI^ = I and 11^11^ = n^TI^ = 0, 
moreover STI^ = ±11$, i.e. the projection TI^ 1 projects onto the ±1 eigenspace ofS. 

Note that if S is involutive, so is —S, which corresponds to the opposite identification of the 
±1 eigenspaces. A matrix K, whose columns are in the span of the +1 eigenspace is said to 
be block-diagonal with respect to the automorphism, while a matrix P, whose columns are 
in the span of the — 1 eigenspace, is said to be 2-cyclic. 

In the context of Lemma[T] we recognize that P £ p is the 2-cyclic part and K 6 t is the 
block-diagonal part. Namely, if X is represented by a matrix, then 

x-( x ~ x ~ + ) 
A ~ \x+- X++ J ' 
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where X' 1 = TI' S XTI J S restricted to the appropriate subspaces. Then, X and X ++ corre- 
sponds to the block-diagonal part K and X h ,X + ~ corresponding to the 2-cyclic part P. 

In passing, we mention that the decomposition l[5]l is called a Cartan decomposition 
whenever the Cartan-Killing form B(X,Y) = tr(adxady) is nondegenerate, hence it can be 
used to introduce a positive bilinear form B^a — —B{X, dc(F)). If this is the case, the linear 
subspaces t,p are orthogonal. 

The involutive automorphism a need not be defined at the group level G and thereafter 
lifted to the algebra by It is possible to proceed the other way around: an involutive 
algebra automorphism da on Q, which automatically produces a decomposition |4|-l|5|, can 
be used to induce a group automorphism a by the relation 

cr(*)=exp(d<T(log*)), (7) 

and a corresponding group factorization Thus, we have an "upstairs-downstairs" view- 
point: the group involutive automorphisms generate corresponding algebra automorphisms 
and vice versa. This "upstairs-downstairs" view is useful: in some cases, the group factor- 
ization Q is difficult to compute starting from x and a, while the splitting at the algebra 
level might be easy to compute from X and da. In other cases, it might be the other way 
around. 



2.2 Lie triple systems 

In Lie group theory Lie algebras are important since they describe infinitesimally the struc- 
ture of the tangent space at the identity. Similarly, Lie triple systems give the structure of the 
tangent space of a symmetric space. 

Definition 2 ([ 14 1) A vector space with a trilinear composition [X, Y,Z] is called a Lie triple 
system (LTS) if the following identities are satisfied: 

(i) [X,X,X] = 0, 

(ii) [X,Y,Z] + [Y,Z,X} + [Z,X,Y]=0, 

(iii) [X, Y, [U, V,W]} = [[X,Y,U], V, W) + [U, [X, Y, V],W] + [U, V, [X,Y,W]\. 

A typical way to construct a LTS is by means of an involutive automorphism of a Lie algebra 
Q. With the same notation as above, the set p is a LTS with the composition 

[X,Y,Z] = [[X,Y],Z]. 

Vice versa, for every LTS there exists a Lie algebra (5 and an involutive automorphism a 
such that the given LTS corresponds to p. The algebra is called standard embedding of 
the LTS. In general, any subset of Q that is closed under the operator 

Tx(-) = ad!(-) = [X,[X,-]] (8) 



is a Lie triple system. It can be shown that being closed under Tx guarantees being closed 
under the triple commutator. 
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3 Application of symmetric spaces in numerical analysis 

3.1 The classical polar decomposition of matrices 

Let Gh(N) be the group of N x N invertible matrices. Consider the map 

a(x)=x- J , xeGL(N). (9) 

It is clear that a is an involutive automorphism of GL(/V). Then, according to Theorem [T] 
the set of symmetric elements G a = {x E GL(N) : o(x) = x~ 1 } is a symmetric space. We 
observe that G a is the set of invertible symmetric matrices. The symmetric space G a is dis- 
connected and particular mention deserves its connected component containing the identity 
matrix /, since it reduces to the set of symmetric positive definite matrices. The subgroup 
G a consists of all orthogonal matrices. The decomposition |2| is the classical polar decom- 
position, any nonsingular matrix can be written as a product of a symmetric matrix and an 
orthogonal matrix. If we restrict the symmetric matrix to the symmetric positive definite 
(spd) matrices, then the decomposition is unique. In standard notation, p is denoted by s 
(spd matrix), while k is denoted by q (orthogonal matrix)]^] At the algebra level, the cor- 
responding splitting is g = p ffi t, where t = {Xe QK N ) '■ da(X) = X} = S0(N), is the 
classical algebra of skew-symmetric matrices, while p = {X 6 Ql{N) : do(X) = —X} is 
the classical set of symmetric matrices. The latter is not a subalgebra of gl(^V) but is closed 
under Tx, hence is a Lie triple system. The decomposition ([SJ is nothing else than the canon- 
ical decomposition of a matrix into its skew- symmetric and symmetric part, X = P + K = 
±(X-do{X)) + t{X+do{X)) = ±(X-X T ) + \(X+X T ).Itis well known that the polar 
decomposition x = sq can be characterized in terms of best approximation properties. The 
orthogonal part q is the best orthogonal approximation of x in any orthogonally invariant 
norm (e.g. 2-norm and Frobenius norm). Other classical polar decompositions x = sq, with 
s Hermitian and q unitary, or with s real and q coinvolutory (i.e. qq = /), can also be fit- 
ted in this framework |2] with the choice of automorphisms a(x) = x~* = x~ T (Hermitian 
adjoint), and a(x) =x respectively (where x denotes the complex conjugate of x). 

The group decomposition x = sq can also be studied via the algebra decomposition. 

3.2 Generalized polar decompositions 

In |21 ] such decompositions are generalized to arbitrary involutive automorphisms, and best 
approximation properties are established for the general case. 

In ll32ll an explicit recurrence is given, if exp(X) = x, exp(S) = s and exp(g) = q then S 
and Q can be expressed in terms of commutators of P and K. The first terms in the expansions 
of S and Q are 

S = P-\[P,K}-\[K,[P,K}) 

+ ±iP,{P,{P,K}}}-±{K,[K,{P,K}}} 

+ [K, [P, [P, [P,K]}}} - ^ [K, [K, [K, [P,K]}}} - ^[[P,K], [P, [P,K]}} + ■■■, (10) 

Q = K-U p APM + Tm l p > i p l p > l p > *]]]] 

+ 720 % IK, [P, [P,K]]}}-m K],[K,[P,K}}}+-- 

2 Usually, matrices are denoted by upper case letters. Here we hold on the convention described in §.2. 
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Clearly, also other types of automorphisms can be considered, generalizing the group 
factorization Q and the algebra splitting {6J. For instance, there is a large source of involu- 
tive automorphisms in the set of involutive inner automorphisms 

o(x) = rxr, da(X) = rXr, (11) 

that can be applied to subgroups of G = GL(n) to obtain a number of interesting factoriza- 
tions. The matrix r has to be involutive, r 2 = I, but it need not be in the group: the factor- 
ization makes sense as long as a(x) is in the group (resp. da(X) is in the algebra). As an 
example, let G = SO(n) be the group of orthogonal matrices, and let r = [— ei,£2, ■ • • ,e n ] = 
I — 2e\e[, where e; denotes the ith canonical unit vector in R". Obviously, r f£ SO(n), as 
detr = — 1; nevertheless, we have (rxr) T (rxr) = r T x T r T rxr = I, as long as x 6 SO(n), since 
r = r T and r T r = I, thus o(x) = rxr 6 SO («). It is straightforward to verify that the subgroup 
G a of Theorem[T]consists of all orthogonal n x n matrices of the form 

_(\ T \ 
q \0q n -i)' 

where q n -\ 6 SO(n— 1). Thus the corresponding symmetric space is G/G a = SO(n) /SO(n — 
1 ) . Matrices belong to the same coset if their first column coincide, thus the symmetric space 
can be identified with the (n — 1 )-sphere S"~ 1 . 

The corresponding splitting of a skew-symmetric matrix V 6 g = 50(n) is 

-(;£H:yms£)-™»«- 

Thus any orthogonal matrix can be expressed as the product of the exponential of a matrix 
in p and one in t. The space p can be identified with the tangent space to the sphere in the 
point (1,0, . . . ,0) T . Different choices of r give different interesting algebra splittings and 
corresponding group factorizations. For instance, by choosing r to be the anti-identity, r = 
[e„,e„_i, . . . ,e\], one obtains an algebra splitting in persymmetric and perskew-symmetric 
matrices. The choice r = [— e\,e%,. . . ,e„, — e„ + i,e„ + 2 . . . ,e2«] for symplectic matrices, gives 
the splitting in lower-dimensional symplectic matrices forming a sub-algebra and a Lie- 
triple system, and so on. In [34, 10], such splittings are used for the efficient approximation 
of the exponential of skew-symmetric, symplectic and zero-trace matrices. In 1121 similar 
ideas are used to construct computationally effective numerical integrators for differential 
equations on Stiefel and Grassman manifolds. 

3.3 Generalized polar coordinates on Lie groups 

A similar framework can be used to obtain coordinates on Lie groups [131, which are of 
interest when solving differential equations of the type 

x = X(t,x)x, jc(0) = e, 

by reducing the problem recursively to spaces of smaller dimension. Recall that x(t) = 
exp(Q (t)) where Q obeys the differential equation Q = dexp^'X, where dexp^ 1 = LJ=o jf ac ^A' 
Bj being the jth Bernoulli number, see J9). 

Decomposing X = TI^X + TI+X e p ffi t, where Tl^X = \(X ± do(X)), the solution 
can be factorized as 

x(t)=exp(P{t))exp(K{t)), 
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where 

P = n-x-[p,n + x} + ^2 2 jc 2jS id 2 P j n-x, <ay=(f|y> (12) 



7=1 



K = dexp^ 1 (n+X - 2 £ {2 2j ~ l)c2jad P J l n~X). (13) 

7=1 

Note that ( p"3} depends on P, however, it is possible to formulate it solely in terms of K, 
n~X and n + X, but the expression becomes less neat. In block form: 

p\ (i o \( ^ -v\fn-z 



where v = adp. The above formula paves the road for a recursive decomposition, by rec- 
ognizing that dexp^ 1 is the dexp -1 function on the restricted sub-algebra t. By introduc- 
ing a sequence of involutive automorphisms (J,-, one induces a sequence of subalgebras, 
g = goDgi D g 2 D . . ., of decreasing dimention, g, + i = Range (iTjr). Note also that the 
functions appearing in the above formulation are all analytic functions of the ad-operator, 
and are either odd or even functions, therefore they can be expressed as functions of ad 2 on 
p. In particular, this means that, as long as we can compute analytic functions of the ad 2 
operator, the above decomposition is computable. 

Thus, the problem is reduced to the computation of analytic functions of the 2-cyclic 
part P as well as analytic functions of adp (trivialized tangent maps and their inverse). The 
following theorem addresses the computation of such functions using the same framework 
of Lemma Q] 

Theorem 2 (| 13 1) Let P be the 2-cyclic part ofX = P+ (X — P) with respect to the involu- 
tion S, i.e. SPS = —P. Let © = P 2 LI^ . For any analytic function \jf(s), we have 

\lf{P) = \lf{0)I+ Wl {&)P + Py x {&)+P W2 (&)P+\lf 2 {&)&, (14) 
where = ^(^(v^) ~ and Vz( s ) = gOKv^) + ¥hV^) ~ 2 V(0))- 

A similar result holds for adp, see 1131 . It is interesting to remark that if 

B T 
A 



then 



Vi(0)/. 

where \)f\ , \f/ 2 as above. In particular, the problem is reduced to computing analytic functions 
of the principal square root of a matrix 1 7 1 : numerical methods that compute these quantities 
accurately and efficiently are very important for competitive numerical algorithms. Typi- 
cally, © is a low-rank matrix, hence computations can be done using eigenvalues and eigen- 
functions of the ad operator restricted to the appropriate space, see 11211 131 . In particular, if 
A, B are vectors, then B T A is a scalar and the formulas become particularly simple. 

These coordinates have interesting applications in control theory. Some early use of 
these generalized Cartan decompositions l(4]i-((5]l (which the author calls 7L 2 -grading) to 
problems with nonholonomic constraints can be found in |1]. In [11], the authors embrace 
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the formalism of symmetric spaces and use orthogonal (Cartan) decompositions with ap- 
plications to NMR spectroscopy and quantum computing, using adjoint orbits as main tool. 
Generally, these decompositions can be found in the literature, but have been applied mostly 
to the cases when [6,6] = {0} or [p,p] = {0}, or both, as 1 12 H 13 I become very simple and 
the infinite sums reduce to one or two terms. The main contribution of 1 13 1 is the derivation 
of such differential equations and the evidence that such equations can be solved efficiently 
using linear algebra tools. See also 1331 for some applications to control theory. 



3.4 Symmetries and reversing symmetries of differential equations 

Let Diff (M) be the group of diffeomorphisms of a manifold M onto itself. We say that a map 
<p e Diff(M) has a symmetry 5? : M — > M if 

j^pj/*- 1 =<p 

(the multiplication indicating the usual composition of maps, i.e. <Pi <p>2 = <Pi ° 92), while if 

M<pMr l = <p _1 , 

we say that Mis a reversing symmetry of 9 [ 1 7 1 . Without further ado, we restrict ourselves to 
involutory symmetries, the main subject of this paper. Symmetries and reversing symmetries 
are very important in the context of dynamical systems and their numerical integration. For 
instance, nongeneric bifurcations can become generic in the presence of symmetries and 
vice versa. Thus, when using the integration time-step as a bifurcation parameter, it is vitally 
important to remain within the smallest possible class of systems. Reversing symmetries, on 
the other hand, give rise to the existence of invariant tori and invariant cylinders 1 19 , 25 28 
[29]. 

It is a classical result that the set of symmetries possess the structure of a group — they 
behave like automorphisms and fixed sets of automorphisms. The group structure, however, 
does not extend to reversing symmetries and fixed points of anti-automorphisms, and in 
the last few years the set of reversing symmetries has received the attention of numerous 
numerical analysts. In 1171 it was observed that the set of fixed points of an involutive anti- 
automorphism srf- was closed under the operation 

<Pi ■ 92 = <Pi 92<Pi € fix -^- > V<pi , <p 2 e fix,c/_ , 

that McLachlan et al. called "sandwich product"[^]lndeed, our initial goal was to understand 
such structures and investigate how they could be used to devise new numerical integrators 
for differential equations with some special geometric properties. We recognise, cfr. j ]2.1| 
that the set of fixed points of an anti-automorphism is a symmetric space. Conversely, any 
connected space of invertible elements closed under the "sandwich product", is the set of 
the fixed points of an involutive automorphism (cfr. Theorem [T} and has associated to it a 
LTS. To show this, consider the well known symmetric BCH formula, 



exp(Z) = exp(X)exp(F)exp(X) 

Z = 2X + Y + \T Y {X) - \T X {Y) + -^X + .- 



3 The authors called the set of vector fields closed under the sandwich product <p t ■ (fo = <pi (fo <Pi a pseu- 
dogroup. 
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1 26 1, which is used extensively in the context of splitting methods 1161 . Because of the 
sandwich-type composition (symmetric space structure), the corresponding Z must be in the 
LTS space, and this explains why it can be written as powers of the double commutator 
operators Tx = ad x ,Ty = ad| applied to X and Y. A natural question to ask is: what is the 
automorphism a having such sandwich-type composition as anti-fixed points? As fix.g/_ = 
{z\o(z) = z }, we see that, by writing z = expZ, we have c(exp(Z)) = (exp(Z)) = 
exp(— Z). In the context of numerical integrators, the automorphism a consists in changing 
the time t to —t. This will be proved in j j3.5| 

If M is a finite dimensional smooth compact manifold, it is well known that the infinite 
dimensional group of Diff(M) of all smooth diffeomorphisms M — > M is a Lie group, with 
Lie algebra Vect(M) of all smooth vector fields on M, with the usual bracket and exponential 
map. It should be noted, however, that the exponential map is not a one-to-one map, not 
even in the neighbourhood the identity element, since there exist diffeomorphisms arbitrary 
close to the identity which are not on any one-parameter subgroup and others which are on 
many. However, the regions where the exponential map is not surjective become smaller and 
smaller the closer we approach the identity 1 24 22 ], and, for our purpose, we can disregard 
these regions and assume that our results are formally true. 

There are two different settings that we can consider in this context. The first is to an- 
alyze the set of differentiable maps that possess a certain symmetry (or a discrete set of 
symmetries). The second is to consider the structure of the set of symmetries of a fixed 
diffeomorphism. The first has a continuous-type structure while the second is more often a 
discrete type symmetric space. 

Proposition 1 The set of diffeomorphisms <p that possess as an (involutive) reversing 
symmetry is a symmetric space of the type G a - 

Proof Denote 

o(<p) =M<pM~ l . 
It is clear that a acts as an automorphism, 

o{<P\9i) = cr((pi)fj(<j()2), 

moreover, if Si is an involution then so is also a. Note that the set of diffeomorphisms <p 
that possess M as a reversing symmetry is the space of symmetric elements G a defined by 
the automorphism a (cfr. Hence the result follows from Theorem[T] 

Proposition 2 The set of reversing symmetries acting on a diffeomorphism <p is a symmetric 
space with the composition S%\ ■ £%2 = Mi. 

Proof If M\ is a symmetry of (p then so is also 1 , since 8%^ (p~ x M\ = (p and the assertion 
follows by taking the inverse on both sides of the equality. In particular, if M\ is a symmetry 
of <p it is also true that is a reversing symmetry of <p~'. Next, we observe that if Si\ 
and S$,2 are two reversing symmetries of <p then so is also S%\,5f~ x S%\, since 

It follows that the composition M\ ■ M% = MiSi^ 1 ^ ^ s an internal operation on the set of 
reversing symmetries of a diffeomorphism (p. 

With the above multiplication, the conditions i)-iii) of Definition [T] are easily verified. 
This proves the assertion in the case when (j) has a discrete set of reversing symmetries. 
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In what follows, we assume that & 6 Diff(M) is differentiable and involutory (X 1 = 
3") and a((p) = &(pX. 

Acting on <p = exp(rX), we have 

doX = — I oexp(tX) = .%X.X, 
dt\t=o 

where .% is the tangent map of & , The pullback is natural with respect to the Jacobi bracket, 

[%XJ,%YJ] = ,%[X,Y]J, 

for all vector fields X, Y. Hence the map da is an involutory algebra automorphism. Let t a 
and p be the eigenspaces of da in Q = diff(M). Then 

= «©P, 

where 

I = {X : %X = X&} 
is the Lie algebra of vector fields that have 5"asa symmetry and 

p = {X : XX = -XX) 

is the Lie triple system, vector fields corresponding to maps that have X as a reversing 
symmetry. Thus, as is the case for matrices, every vector field X can be split into two parts, 

X = X - (x + d<T(X)) + i (X - d<T(X)) = i (x + %X9^ +\{x- %X3^ , 

having X as a symmetry and reversing symmetry respectively. 
In the context of ordinary differential equation, let us consider 

%=F(y), yeR N . (16) 

Given an arbitrary involutive function the vector field F can always be canonically split 
into two components, having M as a symmetry and reversing symmetry respectively. How- 
ever, if one of these components equals zero, then the system l |16| l has 3? as a symmetry or 
a reversing symmetry. 



3.5 Selfadjoint numerical schemes as a symmetric space 



Let us consider the ODE 1 16 1, whose exact flow will be denoted as (p = exp(tF). Backward 
error analysis for ODEs implies that a (consistent) numerical method for the integration of 



16 1 can be interpreted as the sampling at t = h of the flow (pi,(t) of a vector field Ff, (the so 



called modified vector field) which is close to F, 

(Ph(t) = exp(tF h ), Ff, = F + h p E p + h p+1 E p+ i H , 

where p is the order of the method (note that setting t = h, the local truncation error is of 
order h p+1 ). 

Consider next the map o~ on the set of flows depending on the parameter h defined as 



o(<Ph(t)) = <P- h (-t), 



(17) 
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where q>_ h (t) = exp(fF_ A ), with F_ h = F + {-h)PE p + (-h) p+1 E p+1 + ■■■. 

The map o~ is involutive, since a 2 = id, and it is easily verified by means of the BCH 
formula that <j{<P\.h<P2.h) = c( < Pi,/i)o , ( < P2,/i)> hence o~ is an automorphism. Consider next 

Go = {% ■■ e{<Ph) = (Ph 1 }- 

Then <p/, £ G a if and only if <p_/,(— f) = (p^ (t), namely the method (ph is selfadjoint. 

Proposition 3 The set of one-parameter, consistent, selfadjoint numerical schemes is a sym- 
metric space in the sense of Theorem^ generated by O as in \17\ . 

Next, we perform the decomposition We deduce from {T7| that 

da{F h ) = j t \ t=o o( S xp(tF h )) = -(F + (-h)PE p + (-h)"+ l E p+l ) + --- = -F_ h , 
hence, 

t = {F h : da{F h ) = F h } = {F h : -F_ h = F h }, 
is the subalgebra of vector fields that are odd in h, and 

p = {F h : da(F h ) = -F h ) = {F h : =F h }, 

is the LTS of vector fields that possess only even powers of h. Thus, if i 7 /, is the modified 
vector field of a numerical integrator <p/,, its canonical decomposition in t © p is 

F h = \ {F h + da(F h )) + \{F h - da{F h )) 

= \ (F + t (1 - (~l) k )h k E k ) + I (F + £ (1 + (-1)*)^), 

k— p k=p 

the first term containing only odd powers of h and the second only even powers. Then, if the 
numerical method <Ph{h) is selfadjoint, it contains only odd powers of h locally (in perfect 
agreement with classical results on selfadjoint methods l3l). 



3.6 Connections with the generalized Scovel projection for differential equations with 
reversing symmetries 

In 1211 it has been shown that it is possible to generalize the polar decomposition of ma- 
trices to Lie groups endowed with an involutive automorphism: every Lie group element z 
sufficiently close to the identity can be decomposed as z = xy where x 6 G a , the space of 
symmetric elements of a, and y 6 G° , the subgroup of G of elements fixed under <7. Further- 
more, setting z = exp(fZ) and y = exp(7(f)), one has that 7(f) is an odd function of t and 
it is a best approximant to z in G° in G a right-invariant norms constructed by means of the 
Cartan-Killing form, provided that G is semisimple and that the decomposition Q = p © t is 
a Cartan decomposition. 

Assume that <p, the exact flow of the differential equation ( |16| l, has .'M as a reversing 
symmetry (i.e. F 6 p a , where o(<p) = .^(p^ 1 ), while its approximation (p;, has not. We 
perform the polar decomposition 



<P/i = VhXh, o{Vh) = Vh 1 -. = Xh, 



(18) 
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i.e. Xjfh has M as a reversing symmetry, while Xh has M as a symmetry. Since the original 
flow has ^ as a reversing symmetry (and not symmetry), %h is the factor that we wish 
to eliminate. We have = flJ&O^^) -1 . Hence the method obtained composing <p/, with 
o"(<P/?) 1 has the reversing symmetry M every other step. To obtain we need to extract 
the square root of the flow q>i,<y((pi,)~ l . Now, if <j)(t) is a flow, then its square root is simply 
(j)(t/2). However, if (j>i,(t) is the flow of a consistent numerical method (p > 1), namely 
the numerical integrator corresponds to (j>i,(h), it is not possible to evaluate the square root 
(j>i, (h/2) by simple means as it is not the same as the numerical method with half the stepsize, 
^/ 2 (/i/2). The latter, however, offers an approximation to the square root: note that 

h (?) H (?) = ex p(^+^)%) +••• , 

an expansion which, compared with (j>h(h), reveals that the error in approximating the square 
root with the numerical method with half the stepsize is of the order of 

a term that is subsumed in the local truncation error. The choice 

Vh = <Ph/2°(<Ph/2r l = <Ph/2<*(<PHft) (19) 

as an approximation to y/i, (we stress that each flow is now evaluated at t = h/2), yields a 
map that has the reversing symmetry M at each step by design, since 

a(f h ) = a((p h/2 a((p h/2 )) = <7{<p h / 2 )<p h/ l 2 = Wh l - 

Note that iff/, = p^cC^j)' wrlere = ( P*ft/2( — ^ s tne mverse ( or adjoint) 

method of <Ph/2- If ® is given by ( |l7j> , then ^((p^p) = (p^ 2 {h/2) and this algorithm is 
precisely the Scovel projection ll27ll originally proposed to to generate selfadjoint numer- 
ical schemes from an arbitrary integrator, and then generalized to the context of reversing 
symmetries ITTl . 

Proposition 4 The generalized Scovel projection is equivalent to choosing the G a -factor in 
the polar decomposition of a flow (pi, under the involutive automorphism o(<p) = M<pM~ , 
whereby square roots of flows are approximated by numerical methods with half the stepsize. 

3.7 Connection with the Thue-Morse sequence and differential equation with symmetries 

Another algorithm that can be related to the generalized polar decomposition of flows is 
the application of the Thue-Morse sequence to improve the preservation of symmetries by 
means of a numerical integrator [ 8 1 . Given an involutive automorphism a and a numerical 
method <p/, in a group G of numerical integrators, Iserles et al. [8 ] construct the sequence of 
methods 

(pW--<p h , (p^:=(p^a((p [k] ), * = 0,1,2,.... (20) 

Since (pW = a°(p^ k \ it is easily observed that the k-th method corresponds to composing 
CT°<pW and cr'fpW according to the k-th Thue-Morse sequence 01101001 .. ., as displayed 
below in Table [T] (see O0II18I ). Iserles et al. (QD) showed, by a recursive use of the BCH 
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(9)rj0((p) f 7 1 (<p) O - (<p)a ((p) f 7 1 (<p) 
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Table 1 Thue-Morse iterations for the method <pW . 



formula, that each iteration improves by one order the preservation of the symmetry y by 
a consistent numerical method, where 5? is the involutive automorphism such that <j(<j>) = 
y^y^. The main argument of the proof is that if the method has a given symmetry 
error, the symmetry error of a(<j)^) has the opposite sign. Hence the two leading symmetry 
errors cancel and has a symmetry error of one order higher. In other words, if the 

method <p/, preserves J? to order p, then (pl A 'l preserves the symmetry 5? to order p + k every 
2 k steps. As a changes the sign of the symmetry error only, if a method <p/, has a given 
approximation order p, so does cr((p/,). Thus Q>- > = a((p/,) can be used as initial condition in 

l |20[ >, obtaining the conjugate Thue-Morse sequence 10010110 By a similar argument, 

also 

X l0] :=<PhOT X l0] :=a(<Ph), X [k+l] ■= o(Z [ V ] , £ = 0,1,2,..., (21) 
generate sequences with increasing order of preservation of symmetry. 

Example 1 As an illustration of the technique, consider the spatial PDE u t = u x + u y +f(u), 
where / is analytic in u. Let the PDE be defined over the domain [— L,L] x [— L,L] with 
periodic boundary conditions and initial value uq. If uq is symmetric on the domain, i.e. 
uo{y,x) = uo(x,y), so is the solution for all t. The symmetry is a(u(x,y)) = u(y,x). Now, 
assume that we solve the equation by the method of alternating directions, where the method 
(p corresponds to solving with respect to the x variable keeping y fixed, while a((p) with 
respect to the y variable keeping x fixed. The symmetry will typically be broken at the 
first step. Nevertheless, we can get a much more symmetric solution if the sequence of the 
directions obeys the Thue-Morse sequence (iteration <pW) or the equivalent sequence given 
by iteration x ^' ■ This example is illustrated in Figures 



Tfl2 



3.8 Connections with a Yoshida-type composition and differential equations with 
symmetries 

In a famous paper that appeared in 1990 (J3TJ) Yoshida showed how to construct high or- 
der time-symmetric integrators starting from lower order time-symmetric symplectic ones. 
Yoshida showed that, if <p is a selfadjoint numerical integrator of order 2p, then 

<Pah(at)(pp h (pt)(p ah (at) 

is a selfadjoint numerical method of order 2p + 2 provided that the coefficients a and jS 
satisfy the condition 



2a + P = 1 
2a 2p+l + p 2p+l =0, 
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Fig. 1 "Modes" of the symmetry error max^ \u((t,x,y) — u(t,y,x)\ for the method of alternating directions 
forthePDE U[ — u x -f- Uy -\- 2 x 10 u , on a square, with initial condition uq — a x ■ and periodic boundary 
conditions. The integration is performed using the first order method iij <+ \ =u^ + h(Du^ + (Forward 
Euler), where D is a circulant divided difference discretization matrix of the differential operators d x ,d y 
(second order central differences). The experiments are performed with constant stepsize h = 10~ 2 . The 
order of the overall approximation is the same as the original method (first order only), the only difference is 
the order of the directions, chosen according to the Thue-Morse sequence. The different symbols correspond 

to different sampling rates: every step, every second step, fourth, eight, From top left to bottom right: 

sequences '0', '01', '0110', '01101001', etc. 



whose only real solution is 

In the formalism of this paper, time-symmetric methods correspond to G^-type elements 
with a as in and it is clearly seen that the Yoshida technique can be used in general to 
improve the order of approximation of G CT -type elements. 

A similar procedure can be applied to improve the order of the retention of symme- 
tries and not just reversing symmetries. To be more specific, let 5? be a symmetry of the 
given differential equation, namely = F 5^ , with 5? ^ id, Sf~ x = 5?, 5^* denoting the 



pullback of 5? to Q = Vect(M) (see [ 3.4 1. Here, the involutive automorphism is given by 



o(p h (t) = y<p h (t)y, 



so that 



p = {p-. y*p = -py}, z = {K-.y*K = Ky}. 
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Fig. 2 Global error and symmetry error max Xt y\u((t,x,y) — u(t,y,x)\ versus step size h for the method of 
alternating directions for the PDE u, = u x + u y + 2 x 10~ 3 u 2 , as above. The basic method (f> (method "0") is 
now a symmetric composition with the Heun method: half step in the x direction, full step in the y direction, 
half step in the x direction, all of them performed with the Heun method (improved Euler). It is immedi- 
ately observed that the symmetry error is improved by choosing the directions according to the Thue-Morse 
sequence. The order of the method (global error) remains unchanged (order two). 



Proposition 5 Assume that <Ph(t) is the flow of a self-adjoint numerical method of order 2p, 
q> h (t) = exp{tF h ), F h =F + h 2 i , E 2p + h 2 P +2 E 2p +2 + ■ ■ ■ , where Ej = Pj + Kj, Pjep, Kj e t 
and F has 5? as a symmetry. The composition 

•pF'CO = 9ah(at)o(<p bh {bt))<p ah {at), (23) 

with 

a= ^77^ — — T , b=\ — la, 

2 + 2 1 /( 2 P+ 1 ) 
has symmetry error 2p + 2 at t = h. 

Proof Write lf23| as 

<Ph\t) = exp{atF ah )exp(btda{F bh ))exp(atF ah ). 

Application of the symmetric BCH formula, together with the fact that da acts by changing 
the signs on the p-components only, allows us to write the relation ( |2"3"j l as 

= exp((2a + b)tF + (2(at)(ah) 2p + (bt)(bh) 2p )K 2p (24) 
+ {2at{ah) 2p - bt(bh) 2p )P 2p + ff(th 2p+2 ) + ff{Ph 2p ) +■■■), 

where the ff(th 2p+2 ^j comes from the E 2p +2 term and the ff(t^h 2pS ) from the commutation 
of the F and E 2p terms (recall that no first order commutator appears in the symmetric 
BCH formula). The numerical method is obtained letting t = h. We require 2a + b = 1 for 
consistency, and 2a 2p+1 — b 2p+l = to annihilate the coefficient of P 2p , the lowest order 

p-term. The resulting method pi (f) retains the symmetry .y to order 2p + 2, as the first 
leading symmetry error is a G term. 
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This procedure allows us to gain two extra degrees in the retention of symmetry per 
iteration, provided that the underlying method is selfadjoint, compared with the Thue-Morse 
sequence of [8| that yields one extra degree in symmetry per iteration but does not require 
selfadjointness. As for the Yoshida technique, the composition l |23| l can be iterated k times 
to obtain a time-symmetric method of order 2p that retains symmetry to order 2{p + k). The 
disadvantage of \23) , with respect to the classical Yoshida approach is the fact that the order 
of the method is retained (and does not increase by 2 units as the symmetry error). The main 
advantage is that all the steps are positive, in particular the second step, b, whereas the jS is 
always negative for p > 3 and typically larger than a, requiring step-size restrictions for stiff 
methods. In the limit, when p — > °°, a, b — > 1/3, i.e. the step lengths become equal. Thus, the 
proposed technique for improving symmetry is of particular interest in the context of stiff 
problems. This is illustrated by the following example and Figure [3] 

Example 2 Consider the PDE u, = V 2 u — u(u— l) 2 , defined on the square [—1,1] X [—1,1], 
with a gaussian initial condition, m(0) = e ~ 9jr ~ 9 - v ~, and periodic boundary conditions. The 
problem is semi-discretized on a uniform and isotropic mesh with spacing 8 = 0.1 and 
is reduced to the set of ODEs it = (D xx + D yy )u + f{u) = F(u), where D xx = /®I>2 and 
D yy = £>2 ® I, D2 being the circulant matrix of the standard second order divided differences, 
with stencil i [1,-2,1]. For the time integration, we consider the splitting F = F\ + F2, 
where F\ = D xx + 5/ and F2 = D yy +5/ and the second order self-adjoint method^] 

^=C°C o C°C- (25) 

We display the global error and the symmetry error for time-integration step sizes h = 3ho x 
[1, j, . . . , gj], where the parameter ho is chosen to be the largest step size for which the basic 
method <p in ( |25[ l is stable. The factor 3 comes from the fact that both the Yoshida and our 
symmetrization method l |23| l require 3 sub-steps of the basic method. So, one step of the 
Yoshida and our symmetrising composition can be expected to cost the same as the basic 
method $25\ . 

As we can see in Figureji] the Yoshida technique, with a = 1/(2 — 2 1 / 3 ) and jS = 1 — 2a, 
does the job of increasing the order of accuracy of the method from two to four, and so does 
the symmetry error. However, since a < 1/2, the jS-step is negative, and, as a consequence, 
it is observed that the method fails to co nverg e for the two largest values of the step size. 



Conversely, our symmetrization method (23 1 has a = 1/(2 + 2 1 / 3 ) and b = 1 — 2a, with 



b positive, and the method converges for all the time-integration steps. As expected, the 
order is not improved, but the symmetry order is improved by two units. The symmetry c is 
applied by transposing the matrix representation u, -j of the solution at (i8,j8), before and 
after the intermediate Z?-step. Otherwise, the two implementations are identical. 



4 Conclusions and remarks 

In this paper we have shown that the algebraic structure of Lie triple systems and the factor- 
ization properties of symmetric spaces can be used as a tool to: 1) understand and provide a 
unifying approach to the analysis of a number of different algorithms; and 2) devise new al- 
gorithms with special symmetry /reversing symmetry properties in the context of the numer- 
ical solution of differential equations. In particular, we have seen that symmetries are more 

4 We have simply split the nonlinear term in two equal parts. Surely, it can be treated in many different 
ways, but that is besides the illustrative scope of the example. 
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Fig. 3 Error versus step size for the Yoshida technique (left) and the symmetrization technique |23| (right) 
for the basic self-adjoint method \25\ applied to a semidiscretization of u, = V 2 u — u(u — l) 2 . The order with 
the Yoshida techique increases by two units, but the method does not converge for the two largest values of 
the step size. Our symmetrization technique improves by two units the order of retention of symmetry (but 
not global error). Due to the positive step-sizes, it converges also for the two largest values of the step size. 
See text for details. 



difficult to retain (to date, we are not aware of methods that can retain a generic involutive 
symmetry in a finite number of steps), while the situation is simpler for reversing symme- 
tries, which can be achieved in a finite number of steps using the Scovel composition. So 
far, we have considered the most generic setting where the only allowed operations are non- 
commutative compositions of maps (the map <p, its transformed, cr((p), and their inverses). If 
the underlying space is linear and so is the symmetry, i.e. o{(p\ + <f>2) = cr(<pi ) + a((pi), the 
map (pi, = S£±^i£y obviously satisfies the symmetry a, as <j((Ph) = %■ Because of the lin- 
earity and the vector-space property, we can use the same operation as in the tangent space, 
namely we identify a and da. This is in fact the most common symmetrization procedure 
for linear symmetries in linear spaces. For instance, in the context of the alternating-direction 
examples, a common way to resolve the symmetry issue is to first solve, say, the x and the 
y direction, solve the y and the x direction with the same initial condition, and then average 
the two results. 

In this paper we did not mention the use and the development of similar concepts in a 
more strict linear algebra setting^] Some recent works deal with the further understanding 
of scalar products and structured factorizations, and more general computation of matrix 
functions preserving group structures, see for instance [15,5,6] and references therein. Some 
of these topics are covered by other contributions in the present BIT issue, which we strongly 
encourage the reader to read to get a more complete picture of the topic and its applications. 
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5 It seems that numerical linear algebra authors prefer to work with Jordan algebras (see §2), rather than 
Lie triple systems. We believe that the LTS description is natural in the context of differential equations and 
vector fields because it fits very well with the Lie algebra structure of vector fields. 
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